% SIMPSON'S
% ------------------------------------------------------------
%  Input:  f = function
%         df = jacobian of f
%         t0 = initial time
%         tf = final time
%         y0 = inital value
%         h  = step size 
% ------------------------------------------------------------
%  Output: y = array indexed by time step ranging from t0->tf
% ------------------------------------------------------------
% 

function [y,t] = simpson(f, df, t0, tf, y0, h)
  

N = 1 + ceil ((tf-t0)/h - 1e-10);
t = zeros (1,N);
y = zeros (length(y0),N);

t(1) = t0;
y(:,1) = y0;

%initialize 'y_(n+1)' with Euler
if N >= 2
   fyold = feval(f,y(:,1),t(1));
   [y(:,2),t(2)] = euler(f,t0,t0+h,y0, h, 0);
   fynew = feval(f,y(:,2),t(2));
end


% use Paulo-made newton method to estimate 'y_(n+2)';  


for i = 3 : N
   t(i) = t(i-1) + h;
   v = [t(i);-h/3;-(y(:,i-2)+h*(fyold+4*fynew)/3)];
   fyold = fynew;
   [y(:,i),fynew] = newton('f','df',y(:,i-1),v,f,df);
end
%print out exact solution for test2

%yE = zeros(1,N);
%for i = 1 : N
%  yE(i) = exp(-t(i)*sin(2*t(i)));
%end

%plot(t, y(1,:), 'm--', t, yE, 'g');